Evaluating Block Algorithm Variants in LAPACK

The LAPACK software project currently under development is intended to provide a portable linear algebra library for high performance computers. LAPACK will make use of the Level 1, 2, and 3 BLAS to carry out basic operations. A principal focus of this project is to implement blocked versions of a number of algorithms to take advantage of the greater parallelism and improved data locality of the Level 3 BLAS. In this paper, we describe our work with variants of some of these algorithms and the performance data we have collected. LAPACK is planned to be a collection of Fortran 77 subroutines for the analysis and solution of systems of simultaneous linear algebraic equations, linear least-squares problems, and matrix eigenvalue problems [1]. This project will combine the functionality of LINPACK and EISPACK in a single package, incorporate recent algorithmic improvements, and restructure the algorithms to use the Level 2 and 3 BLAS (Basic Linear Algebra Subprograms) for efficiency on today’s high-performance computers. We are investigating variant versions of many of the routines in LAPACK. The building blocks of the LAPACK library are the BLAS, a set of standard subroutines for the most common operations in linear algebra [2,3,4]. The original set of BLAS, consisting of vector-vector operations, was used in LINPACK. Recently, specifications have been drawn up for matrix-vector operations (Level 2 BLAS) and matrix-matrix operations (Level 3 BLAS) to meet the demands of multiprocessing, vectorization, and hierarchical memory in today’s high-performance computers. In particular, the Level 3 BLAS perform O (n 3 ) operations on O (n 2 ) data elements, which helps to improve the ratio of computation to memory references on machines that have a memory hierarchy. This paper describes some of the block factorization routines in LAPACK. The blocked version calls the Level 3 BLAS and, if necessary, an unblocked version of the algorithm to do the processing within a block. The unblocked version calls only Level 1 and 2 BLAS routines and is called directly from the blocked routine if the user has set the blocksize to 1. The LU decomposition is derived by equating the product of a unit lower triangular matrix L and an upper triangular matrix U to the original matrix A . As an illusH A A A I A 31 A 21 A 11